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Abstract 

Recent improvements in an unstructured -grid method for large-scale aerodynamic design are presented. Pre- 
vious work had shown such computations to be prohibitively long in a sequential processing environment. Also, 
robust adjoint solutions and mesh movement procedures were difficult to realize, particularly for viscous flows. 
To overcome these limiting factors, a set of design codes based on a discrete adjoint method is extended to a mul- 
tiprocessor environment using a shared memory approach. A nearly linear speedup is demonstrated, and the 
consistency of the linearizations is shown to remain valid. The full linearization of the residual is used to precon- 
dition the adjoint system, and a significantly improved convergence rate is obtained. A new mesh movement algo- 
rithm is implemented and several advantages over an existing technique are presented. Several design cases are 
shown for turbulent flows in two and three dimensions. 


Introduction 

With the advent of modem computer architectures, aerody- 
namic designers have sought to make use of high-fidelity compu- 
tational fluid dynamics (CFD) codes in their everyday design ef- 
forts. While considerable progress has been made towards this 
goal, realistic use of such tools remains hindered by the extreme 
computational burden associated with such an endeavor. 

A large focus has recently been placed on design algorithms. 
In the area of gradient-based optimization, research has focused 
on several methods for obtaining sensitivity information, and 
many of these approaches rely on an adjoint- van able formulation 
for efficiently computing sensitivity derivatives. The adjoint 
technique is particularly attractive for aerodynamic design prob- 
lems in which there are a large number of design variables, yet 
relatively few constraints. Examples of both continuous and dis- 
crete approaches to this method can be found in Refs. 1-10. 

In Refs. 1-4, a discrete adjoint technique has been imple- 
mented on unstructured grids for two- and three-dimensional 
flows. This work was primarily aimed at performing accurate lin- 
earizations of Reynolds-averaged Navier-Stokes solvers, using 
both compressible and incompressible formulations. Results in- 
dicated highly accurate sensitivity information for fully turbulent 
flows. However, the cost of such computations in a sequential- 
processing environment prevented large-scale design cases from 
being pursued. The preconditioning strategy used for the adjoint 
system in these references was based on a first-order lineariza- 
tion of the residual and often led to poor convergence rates. In 
addition, experience showed that the combination of a distance 
function approach and tension-spring analogy used for mesh 
movement was insufficient when large changes in the geometry 
were necessary. This procedure was also intolerant of initial 
meshes with poor quality. 

In the current work, the linearizations developed in Refs. 1-4 
are modified to run in a parallel processing environment. The do- 
main decomposition and parallelization strategies are discussed, 
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resulting speedups are demonstrated, and the linearizations are 
shown to remain consistent. A new preconditioning strategy for the 
adjoint solver is implemented and significantly improved conver- 
gence is demonstrated for turbulent flow. A new mesh movement 
strategy based on modified linear elasticity theory is also adopted, 
and several advantages over the previous approach are presented. 
Several design cases are also shown. 

Nomenclature 

t> c (i Lift and drag coefficients 

c Chord 

D Vector of design variables 

/ Cost function 

/ Identity matrix 

L Lagrangian function 

Q Vector of dependent variables 

R Discretized residual vector 

t Time 

w, v Nodal displacements 

^ Vector of nodal displacements 

V Volume of control volume 

X Computational mesh 

A Vector of costate variables 

v Poisson’s ratio 

Design Methodology 

Flow Equations 

The governing flow equations are the Reynolds-averaged 
Navier-Stokes equations, 11 coupled with the one-equation turbu- 
lence model of Spalart and Allmaras. 12 The flow solvers used in the 
current work are described at length in Refs. 4, 13, and 14, The 
codes use an implicit, upwind, finite-volume discretization, in 
which the dependent variables are stored at the mesh vertices. In- 
viscid fluxes at cell interfaces are computed using the upwind 
schemes of Roe'\ van Leer 16 , or Osher 17 . Viscous fluxes are 
formed using an approach equivalent to a central -difference Galer- 
kin procedure. Temporal discretization is performed using a back- 
ward-Euler time-stepping scheme. The meshes used in this study 
have been generated using the software described in Refs. 18 and 
19. 

An approximate solution of the linear system of equations 
formed at each time step is obtained using several iterations of a 
point- iterative scheme in which the nodes are updated in an even- 
odd fashion, resulting in a Gauss-Seidel-type method. 



The turbulence model is solved separately from the flow equa- 
tions at each time step, using a backward-Euler time-stepping 
scheme. The resulting linear system is solved using the same point- 
iterative scheme employed for the flow equations. The turbulence 
model is integrated all the way to the wall without the use of wall 
functions. 

Adjoint and Design Equations 

Given a steady-state flow solution in the form ot 
R(D,Q y X) = 0 , a Lagrangian function can be defined as 

L(D, Qy Xy A) = f(Dy Qy X) + A r #?(D, Q , X) 0 ) 

where /(D, g, X) represents a cost function to be minimized and 
A represents a vector of Lagrange multipliers, or costate variables. 
Differentiating this expression yields the following: 




Since the vector of costate variables is essentially arbitrary, the co- 
efficient multiplying [3g/3D] r can be eliminated using the fol- 
lowing equation: 



Eq. 3 represents the discrete adjoint equation for the design prob- 
lem. Once the solution for A has been formed, the remaining terms 
in Eq. 2 can be evaluated to give the desired sensitivity information: 



The adjoint equation given in Eq. 3 represents a linear set of 
equations for the costate variables A . Although this system can be 
solved directly using GMRES, 2 " a time-like derivative is added and 
the solution is obtained by marching in time, much like the flow 
solver: 


where 
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A n+I = A" + A" A 


(5) 


(6) 


The time term can be used to increase the diagonal dominance for 
cases in which GMRES alone would tend to stall. This ultimately 
results in a more robust adjoint solver. 

In Refs. 1-5, an incomplete LU decomposition of the matrix ob- 
tained from a first-order accurate discretization is used to precondi- 
tion the linear system. The preconditioning is applied on the left 
and no fill-in is allowed (ILU[0]). 21 

Domain Decomposition Methodology 

In the current work, the mesh partitioner MeTiS"~ is used to di- 
vide the original mesh into subdomains suitable for a parallel envi- 
ronment. Given the connectivities associated with each node in the 
mesh and the number of partitions desired, MeTiS returns an array 
that designates a partition number for each node in the mesh. The 
user is then responsible for extracting the data structures required 
by the specific application. 


Due to the gradient terms used in the reconstruction procedure, 
achieving second-order accuracy in the flow solver requires infor- 
mation from the neighbors of each mesh point as well as the points 
adjacent to these neighbors. In the present implementation, the gra- 
dients of the dependent variables are first computed on each mesh 
partition and then the results are scattered onto neighboring parti- 
tions. This approach dictates that a single level of “ghost nodes be 
stored on each processor. These ghost nodes that are connected to 
mesh points on the current partition are referred to as “level- T 
nodes. Similarly, the neighbors of level-1 nodes that do not lie on 
the current partition are designated “level-2 nodes. This terminol- 
ogy is illustrated graphically in Fig. 1 . 

The adjoint solver requires similar information: however, unlike 
the flow solver, residual contributions must be written into off-pro- 
cessor memory locations associated with level-2 mesh points. This 
implies that a second level of ghost information must be retained 
along partition boundaries. 

Software has been developed to extract the required information 
from a pre-existing mesh based on the partitioning array provided 
by MeTiS. This domain decomposition operation is done prior to 
performing any computations. The user is also able to read in exist- 
ing subdomains and their corresponding solution files and reparti- 
tion as necessary. This capability is useful in the event that addi- 
tional processors become available or processors currently being 
employed must be surrendered to other users. In addition, software 
has been developed that reassembles partition information into glo- 
bal files and aids in post-processing the solutions. 

Parallelization Strategy 

Each of the codes has been modified to run in a multiprocessor 
environment using a shared memory implementation. This ap- 
proach has been chosen because the primary hardware to be utilized 
is a Silicon Graphics Origin 2000 system. In the current implemen- 
tation, ghost information is exchanged across partition boundaries 
by loading data into global shared arrays which are accessible from 
each processor. Simple compiler directives specific to the Origin 
2000 system are used to spawn child processes for each partition in 
the mesh. This approach scales well and is readily extendable to a 
message-passing or OpenMP 2 ^ implementation. The convergence 
rate of the flow solution is independent of the number of proces- 
sors, whereas the convergence of the adjoint solver varies slightly 
since the preconditioner is only applied locally on each mesh parti- 
tion. 




Number of Processors 


Figure 2. Parallel speedup obtained for the flow solver. 
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Figure 3. Parallel speedup obtained for the adjoint solver. 



Figure 4. Surface mesh for viscous ONERA M6 wing. 
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Figure 5. Location of design variables for ONERA M6 wing. 

Table 1. Sensitivity derivatives for turbulent flow over ONERA 
M6 wing computed in parallel. 


Design 

Variable 

Finite 

Difference 

Adjoint 

Percent Error 

Camber #3 

2.7762 

2.7763 

0 . 004 % 

Thickness #4 

-0.03970 

-0.03971 

0.025% 

Twist #4 

0.00747 

0.00747 

0.000% 

Shear # 1 

0.62023 

n 
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The speedup obtained by parallelizing the flow and adjoint solv- 
ers is demonstrated in Figs. 2 and 3. It can be seen that a nearly lin- 
ear speedup is obtained. For this test, turbulent flow over the 
ONERA M6 wing shown in Fig. 4 is computed. The mesh contains 
359.536 nodes with a wall spacing of 2 x K) -6 of the mean aerody- 
namic chord (MAC). The surface mesh consists of 9.129 nodes. 

To verify that the linearizations performed in Refs. 1-4 have re- 
mained consistent through the port to the parallel environment, sen- 
sitivity derivatives obtained using the parallel solvers on eight pro- 
cessors are compared with centered finite differences. Here, 
turbulent flow over an ONERA M6 wing is computed using a 
freestream Mach number of 0.3, an angle of attack of 2° , and a 
Reynolds number of 5x10* based on the MAC. The mesh used for 
this case consists of 16,391 nodes. All results have been converged 
to machine accuracy, and a step size of 1 x 10 5 has been used for 
the finite-difference computations. For this case, the cost function 
is a linear combination of lift and drag, and the design variables 
generated using the software described in Ref. 24 are depicted in 
Fig. 5. It can be seen from Table 1 that the derivatives are highly 
consistent. 

Adjoint Preconditioning Scheme 

In Rets. 1-5, a preconditioned GMRES algorithm has been used 
to solve Eq. 3. In these references, an incomplete LU-factorization 
with no fill-in allowed IILU(O)] is employed as the preconditioner. 
The factorization is based on the first-order linearization of the re- 
sidual, thereby avoiding excessive storage penalties associated with 
the higher-order stencil for the inviscid fluxes. It has been shown in 
Ref. 4 that the GMRES algorithm may stall and a converged adjoint 
solution may be difficult to obtain using this preconditioner, partic- 
ularly for viscous flows. This has been found to be the case for both 
two- and three-dimensional problems. 
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Iteration 

Figure 6. Convergence of the adjoint solution for different 

preconditioners. 

In an effort to develop a more robust adjoint solver, an improved 
ILU(O) preconditioning technique based on the complete lineariza- 
tion of the residual is employed in the current work. As shown in 
Ref. 25, the additional memory required for storing the complete 
linearization is roughly four times that of the first-order matrix for 
three-dimensional problems. This requirement can be somewhat al- 
leviated by utilizing half-precision storage for these terms. As de- 
scribed in Refs. 3 and 4, the linearizations required for the matrix- 
vector products in the GMRES algorithm are stored for the nearest- 
neighbor terms; these linearizations are also stored in half-precision 
in the current work. Experiments have shown that this strategy 
yields a total memory requirement of about 50% more than the pre- 
vious version of the adjoint solver. 

To demonstrate the improved performance using the higher- 
order preconditioner, adjoint solutions are computed in parallel for 
turbulent flow over the ONERA M6 wing shown in Fig. 4 using 
eight processors. The freestream Mach number is 0.84 , the angle 
of attack is 3.06° , and the Reynolds number is 5x10 based on the 
MAC. For this case. 10 GMRES cycles are used with 10 search di- 
rections and 5 restarts. Results for the first- and second-order pre- 
conditioning strategies are shown in Fig. 6. It can be seen that the 
solver based on the first-order preconditioner fails to converge the 
solution, whereas the method employing the complete linearization 
steadily reduces the residual by nearly five orders ot magnitude. 

Mesh Movement Strategy 

As stated in Refs. 1-5, a combination of a distance function ap- 
proach and a tension-spring analogy has previously been employed 
as a means for modifying volume meshes as the geometric shape is 
changed throughout the design process. It has been found that this 
algorithm lacks the robustness necessary for the design environ- 
ment, particularly for large surface deformations, meshes with 
highly distorted cells, and essentially all three-dimensional geome- 
tries. For this reason, a new approach based on modi tied linear elas- 
ticity theory has been implemented. 

In the approach taken in the current work, it is assumed that the 
computational mesh obeys the isotropic linear elasticity relations 
which take the following form in two dimensions: 26 

v 2 « + * = 0 (7) 

1 - 2vdx 

V 2 v + — I— J-V • * = 0 (8) 

1 - 2 vdy 

where v is Poisson's ratio and the nodal displacement vector is 
given by ^ = ui + vj . Despite the assumption of an isotropic ma- 



Figure 7. Near-field view of baseline mesh. 



Figure 8. Near-field view of mesh after applying distance function/ 
tension-spring analogy. 



Figure 9. Near-field view of mesh after applying modified linear 
elasticity method. 

terial, a spatially-varying value of Poisson s ratio is used in order to 
maintain the physical integrity of highly skewed cells. This value 
has been chosen based on heuristics and is set so that the coefficient 
1/(1 _2v) is equal to the aspect ratio of the local cell, in this 
manner, low aspect ratio cells mimic compressible materials such 
as cork, while high aspect ratio cells tend to behave in an incom- 
pressible fashion, much tike rubber. Since the nodes on the surface 
are constrained, the high aspect ratio cells in the near-wall region 
are not susceptible to compression. A similar mesh movement 
scheme has also been utilized in Ref. 27. Here, anisotropy in Pois- 
son’s ratio is achieved by neglecting the Jacobian associated with 
the transformation between physical and computational coordi- 
nates. In this manner, Poisson's ratio is implicitly determined by 
the cell volumes, so that small cells deform less. 
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Figure 12. Near-field view of mesh with flap rotated using modified 
linear elasticity method. 

To illustrate the advantage over the distance function/tension- 
spnng analogy, the flap on a multielement airfoil has been deflected 
a " d “ ch of the mesh movement strategies has been applied 
Hgure 7 shows a near-field view of the baseline mesh in the region 
between the mam element and flap. Figures 8 and 9 show the 
meshes resulting from the distance function/tension-spring and 
elasticity methods, respectively. It can be seen that the distance 
function/tension-spring analogy allows gaps to form in the mesh. 
voids CaS thC elaStlClty a PP roach Pulls in nearby material to fill the 

It has been found that the elasticity approach also allows for sig- 
nificantly larger geometric deformations. In a similar test, the flap 
ot a multielement airfoil has been deflected from its baseline posi- 
tion shown in Fig. 10. As can be seen from Figs. 1 1 and 12. the dis- 
tance function/tension-spring approach has yielded an invalid 
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Figure 14. Mesh with flap translated Ax/r = 0.02 and rotated 

mesh, while the elasticity formulation has handled the deformation 
in an acceptable manner. Similarly, when a series of flap transla- 
tions and rotations is applied to the geometry shown in Fig. 13 the 
meshes resulting from the elasticity technique maintain a high de- 
gree ot quality as shown in Fig. 14. 

To further quantify the differences between the two mesh move- 
ment schemes, derivatives of lift and drag due to horizontal transla- 
tions of a main element and flap are examined. Ideally, the deriva- 
tive due to a translation of the flap should be equal and opposite in 
sign to a derivative due to an equal and opposite translation of the 
mam element. In practice however, this relationship is affected by 
changes in the topology of the mesh due to the manner in which it 
vanes during a shape modification. 5 

To demonstrate this behavior, derivatives of lift and drag due to 
equal and opposite horizontal translations of the main element and 
flap on the two-element airfoil shown in Fig. 15 have been com- 
puted tor a turbulent flow. For this case, the freestream Mach num- 
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Figure 15. Geometry used for translation derivatives. 


Table 2. Derivatives 
translation using the 


of lift and drag due to flap and main element 
distance function/tension-spring analog}. 


Derivative 

x main 

Xfiap 

M 

* 

.. 

dc/dx 

-1.4785 

2.4033 

0.9248 

dc\/dx 

0.0183 

0.0277 

0.0460 

Table 3. Derh 
translation usi 

natives of lift and drag due to flap ai 
ing linear elasticity. 

[id main element 

Derivative 

*main 

^ flap 

S' 

3c i / 3.x 

-3.8064 

3.8671 

0.0607 

3c, /3x 

0.1722 

-0.1615 

0.0107 


ber is 0.25, the Reynolds number is 9xl0 6 , and the angle of attack 
Table 2 shows the lift and drag derivatives due to translations ol 

the main element and flap using the distance funct.on/tension- 

spring analogy, and the last column is the sum ot these two deriva- 
tives thich ideally would be zero. However, it can be seen from he 
Lie that the derivatives are not at all equal in . and the 

drag derivatives are not even ol opposite sign. This inconsistency 
would be expected to have an adverse effect on an optimization 
nrocedure Table 3 shows the same derivatives obtained using the 
elasticity formulation. Although these derivatives do not sum 
exactlv to zero they do exhibit opposite signs and are much closer 
in magnitude. This tendency has been observed in several cases and 
indicates that the linear elasticity formulation maintains the mesh 
topology in a more consistent fashion. 


Design Cases 


^srr.'sssr--- 

fb^ST. uu*« pressure Uismbunou i, sough, Au esponnren,. 
study of the multielement airfoil geometry shown ,n Fig. 16 has 
been previously performed, and it can be seen from Fig. 17 that 
computations using the baseline geometry are in disagreement™ 
theexperimental results. The model used in the ■ 
non-uniform gap and overlap across the span, and the flap detlMted 
at high dynamic pressures. The goal of the current work is to deter- 
mine a new position of the flap using the pressure distnbution ob- 
tained in the experiment. The improved mesh movement capability 
described abov£ allows for the flap adjustment required by such a 

^TheTreestream Mach number is 0.7 , the angle ot attack is 1.5 , 
and the Reynolds number is 30xl<f . For this case. the design^n- 
ahles are the rotation and x- and y- translations of the p. 

fstgn cycles, the flap has been rotated 3.5“ and reposUioned as 
shown in Fig. 16, although very little change has purred afte he 
first design cycle. It can be seen from Fig. 17 that he result g 
agreement with the experimental results is significantly improved. 
Atoough nol shown. an angle of attack sweep verifies that the new 
position of the flap improves the agreement with the experiment 
across the range of angles of attack. 
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Modified 
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Figure 16. Baseline and modified geometries for multielement 
airfoil problem. 




Figure 18. Density contours for the baseline geometry. 
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Design 
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Figure 20. Cross-sections of the initial and final wing geometries. 


Turbulent Flow Over ONERA M6 Wing 

* ‘ urbU,em fl ° W Wing °P timization performed using the 
ONERA M6 mesh previously shown in Fig. 4. The freestream 
Mach number is 0.84 the angle of attack is 3.06° , and the Rey- 
nolds number is 5x10 based on the MAC. For these conditions 
the baseline geometry exhibits a swept shock extending from the 
root leading edge and a strong normal shock further aft as shown in 
Fig. 18. The objective for this example is to reduce drag while 
maintaining a specified litt. The 20 shape design variables are 
s own in Fig. 19, and the angle of attack is also allowed to vary. 
The design case has been run using approximately 3 days of 
wallclock time on 1 2 processors of an Origin 2000 system. 

Cross-sections of the initial and final geometries can be seen in 
ig. 20. After 5 design cycles, the drag coefficient has been de- 



Design 

Baseline 





Figure 21. Pressure distributions for the initial and final 
geometries. 


wing 



Figure 22. Density contours for the final geometry. 


creased by 15% from 0.0168 to 0.0142. while the lift coefficient has 
maintained its baseline value of 0.253. Pressure distributions at sev- 
eral locations across the span of the wing are shown in Fig. 21 , and 
density contours for the final geometry are shown in Fig. 22. It can 
be seen that the normal shock has been weakened considerably, 
particularly in the outboard section of the wing. 

Turbulent Flow Over Multielement Wing 

In order to handle an arbitrary number of three-dimensional ele- 
ments parameterized by the package described in Ref. 24, software 
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Figure 23. Surface mesh for multielement wing. 
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Figure 25. Cross-sections of the initial and final wing geometries. 


has been developed to combine multiple bodies which employ in- 
dependent parameterizations. To evaluate this capability the base- 
line airfoil depicted in Fig. 16 is extruded in a spanwise direction to 
create a 5° swept wing as shown in Fig. 23. The surface gnd 
shown contains 31,229 nodes and the volume mesh consists ot 
843,385 nodes and 4,796,360 tetrahedra. The adjoint solver re- 
auires roughly 12 gigabytes of storage for the current example. 
q For this 8 caL, the main element and flap are parameterized sepa- 
rately using camber values at the locations shown in Fig. - n a - 
dition to these shape parameters, the deflection as weH as the s verti- 
cal and streamwise positioning of the flap are used as des g 
variables. The angle of attack is also allowed to vary, tor a total of 
34 design variables. The objective is to reduce the drag whilemain- 
taining a specified lift. The freestream Mach number is 0.75, the 
baseline angle of attack is 1 ° , and the Reynolds number is 6.2 mil- 

The design case has been run using 16 processors of an Origin 
2000 system and required approximately 6 days ot wallclock time^ 
After 5 design cycles, the drag coefficient has been reduced trom 
0.0399 to 0.0378, while the lift coefficient has maintained us ongi- 
nal value of 0 437. Cross-sections of the baseline and modified ge- 
ometries can be seen in Fig. 25. while pressure distnbuuons are 
shown in Fig. 26. The flap has been repositioned and can be seen to 
carry a reduced loading, which has been compensated by an in- 
crease in camber across the main element. 



Figure 26. Pressure distributions for the initial and final wing 
geometries. 
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Summary 

An unstructured mesh design methodology based on a discrete 
adjoint formulation has been extended to a multiprocessing envi- 
ronment using domain decomposition and a shared memory ap- 
proach. The parallel implementation has been shown to scale well 
while yielding discretely consistent sensitivity information. 

A preconditioning scheme based on the complete linearization of 
the residual has been demonstrated for adjoint computations. Al- 
though the new strategy requires an increased amount of memory 
due to the larger stencil, it has been found to give superior conver- 
gence rates and hence better reliability. 

An improved mesh movement capability has been developed 
using an approach based on linear elasticity relations. In the current 
work, the scheme is modified to use a spatially-varying value of 
Poisson's ratio to account for highly skewed cells. The new proce- 
dure yields a robust technique which maintains the mesh topology 
in a more consistent fashion than a previous distance function/ten- 
sion-spring analogy. 

Several design examples have been presented which demonstrate 
the improved capability of the current implementation. Reduced 
turnaround time combined with an increased level of robustness has 
enabled previously impractical problems to be addressed. 
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